css2<- read.csv("css05_cv.csv")
css2$part<- css2$cash + css2$vouchers
css2$part[css2$part==2]<-1
css2$part<- css2$part*100


######################################################################################
# Two-Way Fixed Effects Estimators with Heterogeneous Treatment Effects


css2$part<- css2$cash + css2$vouchers
css2$part[css2$part==2]<-1
css2$part<-css2$part*100
css2$voucher_yr<-0
css2$voucher_yr[css2$elec_year>2016]<-1
css2$iswhiteperson<- 0
css2$iswhiteperson[css2$race=="European"]<- 1
css2$iswhiteperson[is.na(css2$race)]<- NA
css2$iswhiteperson[css2$race==""]<- NA
css2$w_vy<- css2$voucher_yr * css2$iswhiteperson
css2$isrep<-0
css2$isrep[css2$party=="Republican"]<- 1
css2$isrep[is.na(css2$party)]<- NA
css2$r_vy<- css2$voucher_yr * css2$isrep
css2$v_vy<- 0
css2$v_vy[css2$sums>12]<- 1
css2$v_vy<- css2$v_vy*css2$voucher_yr 

css2$i_vy<- 0
css2$i_vy[is.na(css2$income)]<- NA
css2$inc2<- as.numeric(substr(css2$income,2,nchar(css2$income)))
css2$incq<- ceiling(ecdf(css2$inc2)(css2$inc2)*4)
css2$i_vy[css2$incq==1 & css2$voucher_yr==1]<- 1
css2$islowinc<-NA
css2$islowinc[css2$incq>1]<-0
css2$islowinc[css2$incq==1]<-1
css2$ishighvot<-NA
css2$ishighvot[css2$sums<13]<-0
css2$ishighvot[css2$sums>12]<-1

#these take a super long time to run on the full sample and sometimes break/force R to abort 
#(at least on my machine)
#so I've commented out a subset to a sample of 50k individuals to optionally run instead
#(50k*9 time periods = 450k observations)
#still does take a few minutes per did_multiplegt call, but more manageable
#results are generally similar but SE can be substantially larger with smaller samples

css3<-css2
set.seed(123)
#ids_to_samp<- unique(css2$id)
#css3<- sample(ids_to_samp, 50000)
#css3<- subset(css2, css2$id %in% css3)
xx1<- did_multiplegt(df= css3, #white
                     Y = "part",
                     G = "id",
                     T = "elec_year",
                     D = "w_vy", brep=4)
xx2<- did_multiplegt(df= css3, #republican
                     Y = "part",
                     G = "id",
                     T = "elec_year",
                     D = "r_vy", brep=4)
xx3<- did_multiplegt(df= css3, #voting
                     Y = "part",
                     G = "id",
                     T = "elec_year",
                     D = "v_vy", brep=4)
xx4<- did_multiplegt(df= css3, #income
                     Y = "part",
                     G = "id",
                     T = "elec_year",
                     D = "i_vy", brep=4)
set.seed(321)
#placebos
xx1_p<- did_multiplegt(df= css3,
                       Y = "part",
                       G = "id",
                       T = "elec_year",
                       D = "w_vy", brep=4, placebo = 1)
set.seed(321)
xx2_p<- did_multiplegt(df= css3,
                       Y = "part",
                       G = "id",
                       T = "elec_year",
                       D = "r_vy", brep=4, placebo = 1)
set.seed(321)
xx3_p<- did_multiplegt(df= css3,
                       Y = "part",
                       G = "id",
                       T = "elec_year",
                       D = "v_vy", brep=4, placebo = 1)
set.seed(321)
xx4_p<- did_multiplegt(df= css3,
                       Y = "part",
                       G = "id",
                       T = "elec_year",
                       D = "i_vy", brep=4, placebo = 1)

# get coefs for left table A26 with: 
xx1$effect
xx1$se_effect
xx2$effect
xx2$se_effect
xx3$effect
xx3$se_effect
xx4$effect
xx4$se_effect

xx1_p$placebo_1
xx1_p$se_placebo_1
xx2_p$placebo_1
xx2_p$se_placebo_1
xx3_p$placebo_1
xx3_p$se_placebo_1
xx4_p$placebo_1
xx4_p$se_placebo_1


############################
# DRDiD 
#this runs pretty quickly

#2015-2019
set.seed(1)
css3<- na.omit(css2)
temptab<- data.frame(table(css3$id))
temptab<- subset(temptab, temptab$Freq==9)
css3<- subset(css3, css3$id %in% temptab$Var1)
y1<- subset(css3$part, css3$elec_year==2019)
y0<- subset(css3$part, css3$elec_year==2015)
D<- subset(css3$w_vy, css3$elec_year==2019)

covariates<- as.matrix(css3[1:NROW(y0),c("isrep","ishighvot","islowinc")])
tt1<- drdid_imp_panel(y1, y0, D, covariates)

D<- subset(css3$r_vy, css3$elec_year==2019)
covariates<- as.matrix(css3[1:NROW(y0),c("iswhiteperson","ishighvot","islowinc")])
tt2<- drdid_imp_panel(y1, y0, D, covariates)

D<- subset(css3$v_vy, css3$elec_year==2019)
covariates<- as.matrix(css3[1:NROW(y0),c("iswhiteperson","isrep","islowinc")])
tt3<- drdid_imp_panel(y1, y0, D, covariates)

D<- subset(css3$i_vy, css3$elec_year==2019)
covariates<- as.matrix(css3[1:NROW(y0),c("iswhiteperson","ishighvot","isrep")])
tt4<- drdid_imp_panel(y1, y0, D, covariates)

tt1
tt2
tt3
tt4

#2013-2017
css3<- na.omit(css2)
temptab<- data.frame(table(css3$id))
temptab<- subset(temptab, temptab$Freq==9)
css3<- subset(css3, css3$id %in% temptab$Var1)
y1<- subset(css3$part, css3$elec_year==2017)
y0<- subset(css3$part, css3$elec_year==2013)
D<- subset(css3$w_vy, css3$elec_year==2017)
covariates<- as.matrix(css3[1:NROW(y0),c("isrep","ishighvot","islowinc")])


tt1<- drdid_imp_panel(y1, y0, D, covariates)
D<- subset(css3$r_vy, css3$elec_year==2017)
covariates<- as.matrix(css3[1:NROW(y0),c("iswhiteperson","ishighvot","islowinc")])
tt2<- drdid_imp_panel(y1, y0, D, covariates)
D<- subset(css3$v_vy, css3$elec_year==2017)
covariates<- as.matrix(css3[1:NROW(y0),c("iswhiteperson","isrep","islowinc")])
tt3<- drdid_imp_panel(y1, y0, D, covariates)
D<- subset(css3$i_vy, css3$elec_year==2017)
covariates<- as.matrix(css3[1:NROW(y0),c("iswhiteperson","ishighvot","isrep")])
tt4<- drdid_imp_panel(y1, y0, D, covariates)
tt1
tt2
tt3
tt4

